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Abstract. Bayesian statistical methods offer a simple and consistent framework for 
incorporating uncertainties into a multi-parameter inference problem. In this work 
we apply these methods to a selection of current direct dark matter searches. We 
consider the simplest scenario of spin-independent elastic WIMP scattering, and infer 
the WIMP mass and cross-section from the experimental data with the essential 
systematic uncertainties folded into the analysis. We find that when uncertainties in 
the scintillation efficiency of XenonfOO have been accounted for, the resulting exclusion 
limit is not sufficiently constraining to rule out the CoGeNT preferred parameter 
region, contrary to previous claims. In the same vein, we also investigate the impact of 
astrophysical uncertainties on the preferred WIMP parameters. We find that within the 
class of smooth and isotropic WIMP velocity distributions, it is difficult to reconcile the 
DAMA and the CoGeNT preferred regions by tweaking the astrophysics parameters 
alone. If we demand compatibility between these experiments, then the inference 
process naturally concludes that a high value for the sodium quenching factor for 
DAMA is preferred. 
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1. Introduction 

Recent years have seen a fervent activity in the direct search of Weakly Interacting 
Massive Particles (WIMPs) in the Galactic dark matter (DM) halo. Besides the well- 
established results of DAMA/Nal and DAMA/Libra [1], which have observed altogether 
13 successive cycles of annual modulation in the nuclear recoil rate consistent with 
the signature of Galactic WIMP scattering, the CoGeNT experiment also claims an 
excess of events that cannot be accounted for by known background sources [2]. If 
interpreted as DM signals, then these results point to a particle mass of few GeVs in 
model independent analyses (e.g., [3-7]), as well as in the frameworks of scalar DM 
(e.g., [8,9]), supersymmetric models (e.g., [10-14]), and hidden sectors (e.g., [15,16]). 

Concurrently, the null results of several other direct detection experiments have led 
to exclusion limits in the WIMP parameter space. For spin-independent scattering, 
CDMS [17,18], XenonlOO [19], XenonlO [20], Edelweiss [21], the CRESST run on 
Tungsten [22], and Zeplin-III [23] have all set relevant bounds. Most notably, the 
limits set by XenonlOO on the WIMP mass and cross-section appear to be incompatible 
with the regions preferred by the DM interpretation of the DAMA and CoGeNT 
results. Prudently though, we note that, depending on the detection techniques, 
direct WIMP searches can be subject to large systematic effects. Indeed, in the case 
of XenonlO, different choices of the scintillation efficiency L e ff can either enhance or 
reduce the compatibility between its exclusion limits and the DAMA/CoGeNT preferred 
parameters [8,24,25]. 

The first goal of this work, therefore, is to address the issue of how to account for 
systematic uncertainties in direct detection experiments. To this end, we employ the 
techniques of Bayesian inference. Bayesian methods provide a simple and consistent 
framework for dealing with nuisance parameters — in this instance, poorly known 
experimental parameters such as L c g — in an inference problem. Once a likelihood 
function has been defined for an experimental result, the nuisance parameters can be 
systematically integrated out of the problem in a procedure known as marginalisation, 
yielding a final posterior probability density function (pdf) for the WIMP parameters 
that incorporates all relevant sources of uncertainties. Because the process of 
marginalisation requires the evaluation of a multi-dimensional integral, Markov Chain 
Monte Carlo (MCMC) methods are particularly well-suited to the purpose. Lastly, 
we note that Bayesian inference is widely used for parameter estimation in precision 
cosmology (e.g., [26]), and has recently also found application in high energy physics, 
e.g., for the exploration of supersymmetric parameter space (e.g., [27-29]), or in view 
of forecasting model expectations for direct DM searches [30-32]. 

In the same vein, the second goal of this work is to incorporate also into the 
picture some degree of uncertainty in the astrophysics. The WIMP-nucleus scattering 
rate in a direct DM search depends on the (unknown) velocity distribution of the 
DM particles in the Galactic halo. Because of its simplicity a common practice is 
to assume a Maxwellian distribution, in which the local DM density, the circular and 
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the escape velocities are fixed at some "standard" values. However, these quantities are 
far from well-constrained by astrophysical observations. Of even less certainty is our 
knowledge of the functional form of the WIMP velocity distribution. Indeed, simulations 
of structure formation suggest that substantial deviations from the isotropic Maxwellian 
form are highly probable (e.g., [33-36]). In this work we investigate several alternative 
velocity distributions. For simplicity we consider only isotropic equilibrium distributions 
consistent with selected spherically symmetric, smooth parametric DM halo density 
profiles motivated by iV-body simulations. Nevertheless, we see no obvious obstacle 
to generalising the analysis also to anisotropic velocity distributions (e.g., [37-42]) and 
non-smooth density profiles (e.g., [43-47]). 

The rest of the paper is organised as follows: after reviewing the basics of direct 
DM searches in section 2, we describe in section 3 various halo profiles and their 
corresponding WIMP velocity distributions. In section 4 we construct the likelihood 
function for each experiment and discuss the modeling of their associated systematics. 
Section 5 contains a detailed explanation of the Bayesian inference procedure. We 
present our inference results in section 6, and conclude in section 7. 



2. The WIMP signal in direct detection experiments 

Direct detection experiments aim to detect or set limits on nuclear recoils arising from 
the scattering of WIMPs off target nuclei. The differential spectrum for such recoils, in 
units of events per time per detector mass per energy, has the form 

dR " e ' dV^| „</(>?(<)), (2.1) 



dE m DM Jv'>v' min dE 
where E is the energy transferred during the collision, p G = Pdm(-Rq) the WIMP density 
in the solar neighbourhood, moM the WIMP mass, da/dE the differential cross section 
for the scattering, and f(v'(t)) is the WIMP velocity distribution in the Earth's rest 
frame normalised such that / d 3 v'f(v'(t)) = 1. The integration in the differential rate 
is performed over all incident particles capable of depositing a recoil energy of E. For 



elastic scattering, this implies a lower integration limit of v' min = ^jM^E /2/a, where 
is the mass of the target nucleus, and \x = tuduMx/ (wdm + M^) is the WIMP-nucleus 
reduced mass. We defer the discussion of the normalised velocity distribution f(v'(t)) 
to section 3. 

The differential cross-section da/dE encodes all the particle and nuclear physics 
information. For coherent elastic scattering it is parameterised as 

da _ M N af (f p Z + (A -Z)f n ) 2 
dE ~ 2fiv» fi 

where \x n = mDM^n/ (wdm + iri n ) is the WIMP-nucleon reduced mass, a^ 1 the spin- 
independent (SI) zero-momentum WIMP-nucleon cross-section, Z (A) the atomic 
(mass) number of the target nucleus used, and f p (f n ) is the WIMP effective coherent 
coupling to the proton (neutron) . We assume the WIMP couples equally to the neutron 



r 2 W, (2.2) 
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and the proton, so that the differential cross-section da/dE is sensitive only to A 2 . The 
nuclear form factor J-(E) characterises the loss of coherence for nonzero momentum 
transfer, and in our analysis we use the Helm form factor [48,49], 

F{E) = 3e -^/2 ^) -fcrcos(AT) > (2 3) 

where s = 1 fm, r = y/R 2 - 5s 2 , R = 1.2 A 1/3 fm, and k = ^/2M N E. 

The total number of recoils expected in a detector of mass M& et in a given observed 
energy range [^1,^2] over an exposure time T is obtained by integrating equation (2.1) 
over energy, 

fall dR 
S(t) = M dct T / dE e(qE) — , (2.4) 

JSi/q dE 

where we have folded into the integral an energy-dependent function e(qE) describing 
the efficiency of the detector. The quenching factor q, defined via S = qE, denotes 
the fraction of recoil energy that is ultimately observed in a specific detection channel 
(scintillation or phonons/heat), and is a detector-dependent quantity. To distinguish S 
from the actual nuclear recoil energy E, the former is usually given in units of keVee 
(electron equivalent keV), while the latter in keVnr (nuclear recoil keV). 



3. The WIMP velocity distribution 

3.1. Halo profiles 

Two astrophysical factors enter into the differential recoil rate (2.1): the local WIMP 
density p & and the corresponding normalised velocity distribution f(v'(t)) in the Earth's 
rest frame (primed v'). These quantities are related via 

p DM (r) = Jd 3 vF(v,r), (3.1) 

where Pdm(^*) is the WIMP density at r from the Galactic Centre (GC) such that 

— * — * 

p Q = pbm(Rq) with R Q = |i? | = 8.5 kpc, and F(v, r) is the WIMP velocity distribution 
in the Galactic frame (unprimed v) whereby f(v'(t)) = F(v, R & )/p & . 

Most analyses in the literature assume a spherically symmetric and isothermal 
distribution for the WIMP around the GC. The WIMP velocities follow the Maxwellian 
distribution F(v,r) ~ exp(— v 2 /v 2 ), where v = 220 km s _1 is the mean velocity in the 
Galactic frame, and the distribution is cut off at v > v csc = 544 km s" 1 . The resulting 
density profile scales as r~ 2 [50], and is normalised to p Q = 0.3 GeV cm" 3 . This is 
known as the Standard Model Halo (SMH). 

However, isothermal DM density profiles are rarely if ever encountered in iV-body 
simulations. Indeed, most simulations find dark matter halos that are often "cuspy" 
and have density profiles that fall off faster than the r~ 2 dependence of the SMH at 
large r. In view of the uncertainty in the exact DM distribution, we consider, besides 
the SMH, four other spherically symmetric DM density profiles found in the literature: 
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(i) Cored isothermal A variant of the SMH, this density profile has the form 



Pdm(^) = Ps 



1 + 



(3.2) 



Unlike the SMH in which the density Pdm(^) diverges as r — > 0, the cored isothermal 
halo has a finite density core whose size and density are characterised by the 
parameters r s and p s respectively. The profile's large r (i.e., r ^> r s ) behaviour, 
however, is similarly to that of the SMH. 

(ii) Navarro-Frenk-White (NFW) Based on iV-body simulation results, Navarro, 
Frenk and White suggested as a universal form for the DM density profile across a 
wide range of halo masses (10 11 — > 1O 15 M ) [51], 



ty* \ 1 / Z' 7" ^ ^ ^ 



A»(r)=^-j ■ (3-3) 

The density here falls off as r -3 at r ^> r s , while at r <C r s we find a r _1 behaviour 
(i.e., the cusp). The NFW profile is formally divergent as r — > 0. For numerical 
stability, however, we introduce in the profile a small core of size e<r s . A related 
density profile is the Moore profile [52], which also exhibits a r~ 3 behaviour at 
r 3> r S) but has a steeper cusp that scales as r -15 . We do not consider the Moore 
profile here because of its similarity to the NFW profile (see section 3.2). 

(iii) Einasto Some recent studies find that the Einasto profile [53] provides as good a 
fit as the NFW profile to DM halos found in iV-body simulations of the concordance 
ACDM model [54]. The Einasto profile has the form 



Pdm(^) = p s exp 



2 17 r 
a IAr f 

where a = 0.17, and its central density is finite 
(iv) Burkert The Burkert profile, 



(3.4) 



pMr)=Ps (l + y (l + ^j , (3-5) 

is a cored profile that appears to provide a good fit to the DM distribution of dwarf 
galaxies [55]. 

All four profiles depend on two parameters p s and r s . However, it is equally valid, 
and perhaps more enlightening, to adopt a parameterisation in terms of the virial mass 
M v j r of the DM halo — defined as the mass contained in a sphere of radius r v j r whose 
average density is 200 times the critical density — and a concentration parameter given 
by c v i r = r vir /r s . The advantage of this parameterisation is that, firstly, it is possible to 
specify directly a prior for M vir based on what we know about the mass of the Milky Way 
from satellite kinematics etc. Secondly, the concentration parameter c vir is well studied 
in iV-body simulations, which again allow us to impose a prior on c v i r in a meaningful 
way. We show how each density profile (3.2) to (3.5) can be expressed in terms of M vir 
and c v i r in Appendix A. 
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3.2. Extracting the velocity distribution 

Given a DM density profile, the underlying DM velocity distribution can be extracted 
by inverting equation (3.1) under the assumption of hydrostatic equilibrium. For 
a spherically symmetric density distribution and assuming an isotropic velocity 
distribution F(v, r) = F(e) in the Galactic frame that depends only on the relative 
energy e = \v 2 > of the system, the solution is given by the Eddington formula [50], 



F(e) 



£ d 2 p DM d^ 1 ( dprjM 



*=o. 



(3.6) 



d^ 2 ^rr^j ^ \ &m 

The function \&(r) is the gravitational potential generated by the DM halo and the 
baryonic matter residing in the Galactic disk and bulge, defined so that \I/(r — > oo) = 0, 
and v esc (r) = ^2^{r) is the escape velocity at r. It is obtained by solving the Poisson 
equation, 

d 2 \I> 2 d^ 

dr 2 " + r dr = ~^ nG ^ PDM + Pdisk + ^ bul § e ]' ( 3 - 7 ) 
where the disk density distribution is given by [30] 

M disk e~ r / rdisk 

Pdisk(r) = - — I , (3.8) 

with M disk = 5 x 10 10 M and r disk = 4 kpc, and the bulge is modelled as a point mass 
sitting at r — 0, 

PbuigcM = M hulge 5 { p\r), (3.9) 

where Mb u i ge = 1-5 x 1O 1O M , and 5p\f) is the 3-dimensional Dirac delta distribution. 

At any given point r, the Eddington formula (3.6) returns a positive and nonzero 
solution for F(e) only up to the escape velocity t> esc at that point. For v > v csc , F(e) is 
by definition zero. Furthermore, the formula shows that the DM velocity distribution 
at R Q , F(ty Q — \v 2 ) 1 depends only on the DM density distribution at r > R Q . Thus, 
halo density profiles sharing the same large r behaviour will yield similar solutions for 
F(^ — \v 2 ) [56]. For this reason the NFW and the Moore profiles are for our purposes 
equivalent (see section 3.1). 

The last step is to rewrite the velocity integral in the differential recoil rate (2.1) 
in terms of F(^ Q — \v 2 ), that is, 

L, dV ^- 2 ^X>,. -5^). (3.10) 



mm 



with 



y 2 = lp , - |2 _ „./2 ..2 



v®\ =v +v !£ + 2v'v S) a, 

v® = |w + w"e,rot| = v Q +u£ ilot cos7cos[27r(t- t )/T] , (3.11) 

where %, and % are, respectively, the Earth's and the sun's velocity in the Galactic 
frame, i>"©, ro t is the Earth's rotational velocity around the sun in the sun's rest frame, 
and 7 = 60° is the inclination of the Earth's rotation plane with respect the the Galactic 
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plane. For our analysis we take v'^ Iot = 29.8 km s 1 , and v = v + 12 km s 1 [49, 57], 
where 



v = \ -r 



dr 



(3.12) 

r=R Q 

is the circular velocity of the local standard of rest. In the time-dependent piece the 
period T is one year, while to corresponds to June 2, the day on which v® reaches its 
maximum. 

Finally, we note that implicit in the recoil rate (2.1) are four astrophysical 
observables: the local DM density p , the Milky Way virial mass M v - ir , the circular 
velocity of the local standard of rest v defined in equation (3.12), and the local escape 
velocity, 

v csc = y^\ r=Re ■ (3.13) 

These can be independently constrained using observations of stellar and satellite 
kinematics. We discuss this point in more detail in section 4.7. 



4. Experiments and their likelihood functions 

The likelihood function C(X\9) denotes the probability of the data X given some 
theoretical prediction 9, and plays a central role in Bayesian inference. In this section 
we describe the likelihood function used for each experiment, as well as the modelling 
of potential systematics. An in-depth discussion of Bayesian methods is deferred to 
section 5. Table 1 summarises the free (MCMC) parameters of our analysis. 



4.1. CD MS St 

The cryogenic CDMS experiment at the Soudan Underground Laboratory operates 
germanium and silicon solid-sate detectors. Two events were observed at 55 and 
95 keVnr in the silicon run (CDMSSi hereafter) on a 0.1 kg detector in an exposure 
of 65.8 kg-days, compatible with an expected background of B n = 3.6 neutrons and 
B e = 0.8 ± 0.6 electrons in the 5 — > 100 keVnr detection window [17]. No quenching 
factor is required for the CDMS experiment, i.e., q — 1. For details of the detector 
efficiency e(qE) we refer the reader to, e.g., [58]. 

We model the corresponding likelihood function with a Poisson(2) distribution,! 

ln£ C DMSSi(2|S,5) = -S-5 + 2 + 21n(^^) , (4.1) 

where S is the expected WIMP signal in the detection window, and B — B n + B e the 
expected background. The likelihood function (4.1) is normalised such that ln£ = if 
the sum of the expected signal and background matches exactly the number of observed 
events. 

| The notation Poisson(n) denotes the Poisson distribution for n observed events. 
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Table 1. MCMC parameters and priors for the WIMP parameter space and 
experimental systematics (nuisance parameters). All priors are uniform over the 
indicated range. 



Experiment 


MCMC parameter 


Prior 


All 


log(m DM /GeV) 


0^3 


All 




_44(_46) -> -38 


DAMA 


<?Na 


0.2 ->■ 0.4 


DAMA 


Qi 


0.06 ->■ 0.1 


XenonlOO 


m 


-0.01 ->■ 0.18 


CoGeNT 


C 


^ 10 cpd/kg/keV 


CoGeNT 


s 


^ 30 keV 


CoGeNT 


G n 


^ 10 cpd/kg/keV 


CDMSGe(LE) 


a 


-0.60 ->■ -0.18 



Since the expected background rate comes with an uncertainty — in this instance, 
B = B ± <jb = 4.4 ± 0.6, it is useful to construct an effective likelihood function C eS by 
marginalising over the background B, 



''-CDMSS 



where 



p(B) 



i (2\S)= / dB C C DMSSi(2\S, B)p(B), 
Jo 

(B — B) 2 ' 



exp 



2al 



(4.2) 



(4.3) 



is the probability density function of B (modelled as a Gaussian distribution). In the 
small cr B /B limit, the resulting effective likelihood has the form, 



m ^ci)MSSi — —S 



B + -# + 2 + In 



a% + (S + B - alf 



(4.4) 



which we use in our inference analysis. 



4.2. CD MSG e 

For their germanium run, the CDMS-II (CDMSGe hereafter) reported two events at 12.3 
and 15.5 keVnr in the 10 — > 100 keVnr window in a total exposure of 612 kg-days [18]. 
The total expected background in the same time frame is B = 0.8 ± 0.1 ± 0.2. For 
our analysis, however, we adopt the fitting formula for the differential background rate 
provided by [24] , 



dN B 
dE 



-0.00295 + 0.463 



/ keVnr\ 



V E 



/(612 kg days) , 



(4.5) 



§ We adopt the small cfb/B limit results whenever B > 3ctb- 
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where the rate has been normalised to B — 0.8 events over the detection window in an 
exposure of 612 kg-days. Exploiting this spectral information, we model the likelihood 
function as a product of two Poisson(l) distributions (for those energies with one event 
each) and a series of Poisson(O) distributions (for those energies with no events) [59], 
that is, 



ln£cDMSG « = - S - B+2+ E ln g + |^) +c „ 



(4.6) 



where S and B are, respectively, the total expected signal and background in the 
detection window, E lf2 = 12.3, 15.5 keVnr, and C norm = J2i=i 2 hi[Md e tTe(gi?j)] is a 
normalisation factor following from the normalisation of the individual Poisson(l) and 
Poisson(O) distributions. See discussion after equation (4.1). 

Marginalising over the total background B (but not the spectral shape) in the 
manner of equations (4.2) and (4.3), we find in the small <Jb/B limit an effective 
likelihood 



^2 

In £r;r)MSr!p — S B H - h 2 + C norm + 



In 




B-a% dN B \ 2 1 dN E 
B dEi) <JB l il 2 BdE l 



(4.7) 



We include in the analysis also null results from three previous searches with the CDMS 
germanium detector, with exposures of 34 kg-days [17], 19.4 kg-days [60], and 397.8 kg- 
days [61], bringing the total exposure to 1063.2 kg-days. The expected background and 
its uncertainty are scaled correspondingly to B — 1.39 and o~b = 0.38 respectively. We 
model the detector efficiency e(qE) after [24]. 



Low energy CDMS The CDMS collaboration has recently re-analysed their germanium 
data — both on their own, and in combination with data from the silicon detectors — with 
a lower energy threshold [62, 63], thereby increasing the experiment's sensitivity to light 
WIMPs. In reference [62], data from 8 germanium detectors (CDMSGe(LE) hereafter) 
were re-analysed using a threshold of 2 keVnr, compared to 10 keVnr in the standard 
analysis. For each detector, the collaboration provides the event energies and the raw 
exposure. After summing up all contributions and applying the efficiency cuts one finds 
a total of 427 counts for 214 kg-days, distributed in the energy range 2 — > 100 keVnr. 
In our analysis, we bin the data in such a way that 16 bins are contained in the energy 
range 2 ->■ 10 keVnr, and 9 in 10 — >■ 100 keVnr. 

A lower energy threshold, unfortunately, is traded at the cost of an increased 
acceptance of background events, because at these low energies the ability of the 
experiment to discriminate between nuclear and electron recoils degrades and the 
ionization signal becomes dominated by noise. Indeed, while the background due to 
surface events, "zero-charge" events, and leakage events are reasonably well known at 
energies > 5 keVnr, between 2 keVnr and 5 keVnr the CDMS collaboration has to rely 
on an extrapolation to model these events in their analysis, as described in figure 1 
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of [62]. Another potential issue is the calibration of the recoil energy near threshold, 
since the ionisation signal is missing. 

Given these considerations, we model the differential background rate in our 
analysis as 

m (F)-i ™ B(i?) ' £>5keVnr, 

Bl ] ~ \ 0.1 X l ^/ ke Vnr)-5] ; 2 < £ /keVnr < 5> ^ 

where m-B(E) corresponds to the black curve in figure 1 of [62], which, for energies above 
5 keVnr, we regard as reliable and free of systematics. For energies below 5 keVnr we use 
an extrapolation function, but allow the slope a to vary subject to a Gaussian constraint 

lnC mB = - { ^^, (4.9) 

where the "best-fit" a = —0.36 reproduces the black curve in figure 1 of [62], and 
o~ a = 0.2 a, chosen based on the error bars in the same figure at 5 keV. 

The expected signal rate in the ith energy bin is then a sum of the DM signal and 
the background rate, 



1 r E i 
' AEt J Ei - 



Ei+AEi/2 

dE 

AEi/2 



lE +ra B (E) 



(4.10) 



where AEi. The likelihood function is given by 



AW /„ ~obs\2 



g? )^ 

In £ C DMSGc(LE) = - X! 9 2 ^ l n ^m B • (4-H) 

i=l Z<J i 

where s° bs is the observed rate in the ith bin, and Oi is the associated error. 

We do not consider the re-analysis of the combined data on the germanium and 
the silicon towers presented in [63] , because the lack of knowledge about the low-energy 
background makes it difficult for us to model the likelihood function. 

4.3. CoGeNT 

The CoGeNT experiment, an ultra low-noise (and hence low-threshold: 0.4 keVee) 
germanium cryogenic detector running at the Soudan Mine, found in a total exposure of 
18.48 kg-days an excess at low energies that cannot be attributed to known background 
sources [2]. Using the energy binning in figure 3 of [2] in the 0.4 — > 3.2 keVee energy 
range, we model the likelihood function as a sum of Poisson(Xj) distributions, 

' Si + bi + r. 



ln£ CG = 



56 

K It ----- / 

x ■ < 4 ' 12 > 



n + X + X in 



i=l L \ 

where X is the number of events observed in the ith. energy bin, is the expected 
signal computed from equation (2.4) with S\ and £2 corresponding respectively to the 
lower and upper energy limits of the bin concerned, and bi and r\ are two background 
components. 
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For the first component hi, we model the differential rate as an exponentially 
decaying function, 

d ^=Cexp(-S/S ), (4.13) 

where C and S are two free parameters. The second component r\ denotes events due 
to two radiation peaks from 65 Zn and 68 Ge decays, whose differential rates are modelled 
as Gaussians centered on the energies £ Zn = 1.1 keVee and S Ge = 1.29 keVee, with a 
common standard deviation fixed by the energy resolution of the detector AS, i.e., 

^=^,z n exp^-^ A ^j, (4.14) 

and similarly for dN Ge /dE. We fix the ratio of the two peak heights to G nj z n /G n ,Ge — 0.7, 
and vary only G n = G n ,Ge- For details about cosmogenic backgrounds we refer 



to [64]. The energy resolution is given by A£/eVee = + F (S/eVee), with 
o~ n = 69.4 and F = 0.29 [2], while the energy-dependent quenching factor is taken to be 
q = 2/[l + yjl + 15.55 (keVee/£)], following [24]. 

Because the background model parameters C , Sq and G n enter into the likelihood 
function (4.12) in a nontrivial fashion, analytical marginalisation in the manner of 
equation (4.2) is cumbersome if not impossible. We therefore treat these parameters as 
MCMC parameters. 



4-4- Xenon 

Xenon is a two phase (liquid/gas) xenon experiment running at Laboratori Nazionali 
del Gran Sasso (LNGS). A nuclear recoil from particle scattering is inferred from the 
simultaneous measurements of scintillation light and ionisation electrons, together with 
the arrival direction. The amount of nuclear recoil energy going into the primary 
scintillation signal is expressed in terms of the number of photoelectrons (PE) produced 
Si, which is related to the nuclear recoil energy E through the relation 

S 1 (E) = L cS (E)L y E^, (4.15) 

>Jee 

where L e g(E) is the energy-dependent scintillation efficiency, L y = 2.2 PE/keVee the 
scintillation efficiency of nuclear recoils relative to that of the 122 keVee 7-rays at 
zero field, and the quantities S nrtSe = 0.95,0.58 denote respectively the electric field 
scintillation quenching factors for nuclear and electron recoils. 

Between 2006 and 2007, the XenonlO collaboration found 13 events for an expected 
background of 7 events in an exposure of 316.4 kg-days in the 2.0 — > 75.0 keVnr 
window [20,65]. The XenonlOO experiment recently released an analysis of 100.9 live 
days of data acquired in 2010, which found three candidate events for an expected 
background of B = 1.8 ± 0.6 in an exposure of 1481 kg days. [19]. Because of this large 
exposure, we consider in this work only the results of XenonlOO from the aforementioned 
data release, although the analysis techniques can easily be generalised for used with 
XenonlO. 



Bayes and present DM direct detection 



12 



As in the case of CDMSGe we include both the total rate and spectral information 
in the likelihood, 



In £ E vents = - S - B + 3 + ^ In ( ^ 

i=i V ^ 



B dN B 

+ ^ 



B dSi 



+ C norm , (4.16) 



where three events are at 8, 20, and 23 PE, respectively, and C norm = Si=i,2,3 m (M det T) 
is a constant normalisation factor. The background has been shown to consist mainly of 
electron recoils, with a flat distribution in energy over the experimental range according 
to measurements and Monte Carlo simulations [66]. We therefore model the differential 
background as dNs/dSi = 0.069/(1481 kg days), normalised to B — 1.8 in the detection 
window of 4 ->■ 30 PE in 1481 kg days. 

The expected WIMP signal S is computed as follows. Firstly, we note that the 
actual conversion between the nuclear recoil energy and the number of PEs produced 
is not deterministic. Given some recoil energy E, the number of PEs produced Si is 
subject to Poisson fluctuations, with equation (4.15) expressing only the expectation 
value Si. In physical terms this means recoils below the nominal energy threshold have 
a finite probability of leaking into the detection window of the experiment, and this 
effect is important for the detection of light WIMPs. The expected number of WIMP 
events as a function of the (discrete) number of PEs generated can be written as 

dR r°° dR 

' dE ^ x PiSilS^E)) , (4.17) 



dSi Jo dE 

where P(Si\Si(E)) denotes a Poisson(n) distribution with expectation Si(E). Summing 
over all PE counts, the total number of events expected in the detector is 



PE max J D 

S = M det T £ (4.18) 
„_ PE . O-Ol 

where for XenonlOO, PE min = 4 and PE max = 30. 

It remains to specify the energy-dependent scintillation efficiency L e g(E). In this 
work we use, 

' L cS (E), E>3 keVnr, 

k max{m[ln(£/keVnr)-ln3] + 0.09, 0}, 1 < E/keVnr < 3. 

(4.19) 

Here, L e g(E) corresponds to the best-fit in figure 1 of [19], which, at E > 3 keVnr, 
is well-constrained by direct measurements. No direct measurements exists at 1 < 
E/keVm < 3, and the "best-fit" provided by the XenonlOO collaboration is merely 
an extrapolation. We therefore treat the extrapolation slope m as a variable, MCMC 
parameter, subject to a Gaussian constraint of 



L cff (£) 



(m — m) 



2 



ln£ Lcff = ~ v 2 ' , (4.20) 

£0 m 



m 



where m = 0.082 reproduces the "best-fit" of [19] in the 1 < E/keVm < 3 region, and 
a rn = 0.04, chosen so that the 2a region coincides approximately with the light blue 
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band. We further restrict m to lie within the range [—0.01,0.18], so that L e s(E) never 
exceeds 0.1 at E — 1 keVnr, or drops to zero at energies above 2 keVnr. 

We note that a somewhat different parameterisation for the uncertainty in L e g(E) 
was adopted in [66], which also accounts for errors in the direct measurements of L e ff 
at E > 3 keVnr. Nonetheless we choose the parameterisation (4.19) because it highlights 
the role of L e g(E) in the low energy region. Indeed, from equation (4.17), we see that the 
main uncertainty in the expected event rate at energies close to the threshold comes from 
Poisson fluctuations in the number of PEs produced, which in turn depend sensitively 
on the unknown slope of L eS (E) at E < 3 keVnr via equation (4.15). The small errors 
around the best-fit L c g(E) at E > 3 keVnr, on the other hand, has little impact on 
the physics close to the threshold and hence also the exclusion power of the experiment 
for light WIMPs. Note that our likelihood function (4.16) automatically takes into 
account uncertainties in L e q(E) in both the total number of signal events and their 
energy dependence, in contrast to the approach of [66], which explicitly assumes the 
(normalised) signal energy spectrum to be independent of these uncertainties. 

Thus, the full likelihood function describing the XenonlOO experiment is 

In £ X enon = In ^Events + h ^L cB , (4.21) 

which we further marginalise numerically over the background events B (in the In ^Events 
term) as per equations (4.2) and (4.3), yielding an effective likelihood that depends only 
on m DM , a^ 1 , and the systematics nuisance parameter m. 

The XenonlO collaboration recently published a low-energy analysis based on 
the ionisation signal (the S 2 signal) [67]. This alternative approach removes the 
dependence of the result on the uncertainties of the scintillation efficiency, thereby 
lowering the detector threshold significantly down to 1 keV. However, without a 
reliable parameterisation of the estimated background, it is not possible to construct 
a meaningful likelihood function for our Bayesian analysis. Since the XenonlO 
collaboration does not provide the necessary information, we refrain from using their 
data in this work. 



4.5. DAM A 

The DAM A/Libra [1,68] experiment at LNGS uses Nal(Tl) crystal radio-pure 
scintillators as targets. The signature of WIMP interactions consists of an annual 
modulation of the signal due to the motion of the Earth through the Galactic halo as 
discussed in section 3.2. For a cumulative exposure of 1.17 ton-year, the collaboration 
reported a positive detection at 8.9a significance. 

For isotropic WIMP velocity distributions such as those considered in this work, the 
expected modulation signal averaged over an observed energy interval [£i,<?2] is given 
by 

1 f S,/ qx I \dR Xn di? x 
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Table 2. Additional MCMC parameters and (uniform) priors related to the modelling 
of the WIMP velocity distribution. 



Density profile 


MCMC parameter 


Prior 


All 


M vir 


1 ->■ 5 x 10 12 M Q 


NFW, Einasto 


C v ; r 


5 ^ 20 [74] 


Cored isothermal, Burkert 


C v ; r 


50 ->■ 200 [75] 



where wx = M x /(M Nll + Mi), and we have explicitly ignored the small contribution 
from channelling [69]. The likelihood follows a Gaussian distribution, 

AW („ _ ^obs\2 

In £dama = - E o 2 > ( 4 - 23 ) 

where and s° hs are the theoretical and the mean observed modulation respectively in 
the iih energy bin, o~i is the associated uncertainty in the observed signal, and we use in 
this analysis the 36-bin data from figure 9 of [68] . The quenching factors q^ a and q\ are 
taken to be free parameters in our analysis, which we vary, respectively, over a range 
representative of the diverse measured values found in the literature [70-73]. 

4-6. Other experiments 

Even though we will not consider their results in our analysis, let us also mention the 
following direct detection experiments. 

The CRESST collaboration has found 32 events on oxygen given an expected 
background of 8.7 ± 1.4 [76]. If interpreted as a WIMP signal, this would point to 
a low mass WIMP with a coherent cross-section in the ballpark of the regions preferred 
by DAMA and CoGeNT data [7]. 

The Edelweiss collaboration recently published the final analysis for their second 
run, reporting 5 events for an expected background of 3, 4 of which are close to the 
threshold [21]. An exclusion bound was set, which, because of the smaller exposure of 
the experiment, is not competitive with those derived from other experiments considered 
in this work. A recent combined analysis of Edelweiss and CDMS has improved the 
sensitivity: a tighter exclusion bound for DM masses above 200 GeV was found, relative 
to the limits obtained by the individual experiment alone [77]. However, this new limit 
is still less constraining than that derived from the XenonlOO experiment. 

4-7. Astrophysics 

In addition to the WIMP mass, cross-section, and the nuisance parameters of the 
direct search experiments, two further free parameters are used to characterise the 
WIMP velocity distribution: the virial mass of the DM halo, and its concentration (see 
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Table 3. Astrophysical constraints on the DM halo profile and the WIMP velocity 
distribution. 



Observable Constraint 

Local standard of rest v° hs = 230 ± 24.4 km s" 1 [78, 79] 

Escape velocity < bs = 544 ± 39 km s _1 [80, 81] 

Local DM density p° bs = 0.4 ± 0.2 GeV cm" 3 [82, 83] 

Virial mass M* s = 2.7 ± 0.3 x 1O 12 M [84, 85] 



table 2). These additional parameters are, however, also constrained by astrophysical 
observations. For this reason, we define a likelihood function for the astrophysics, 

(vo-v^l (v csc -& (P0-P0 8 ) 2 (M vir - M° bs ) 2 

(4.24) 

where the measured values of the various astrophysical observables (see section 3.2 for 
their definitions) and their uncertainties are given in table 3. Note that none of the 
constraints in table 3 assumes a specific parameterisation of the halo profile, which 
allows us to apply them to all halo models we are considering here without running the 
risk of double-fitting. 



5. Statistical inference 



Having specified a theoretical model with free parameters 9 in sections 2 and 3, and 
defined the likelihood functions C(X\6) in section 4, one final step remains to be taken in 
the analysis of the data X: the inference of the posterior probability density as a function 
of the parameters, V(9\X). The posterior pdf represents our state of knowledge about 
the parameters after taking into account the information contained in the data, and has 
an intuitive and straightforward interpretation in that J v V(9\X)d9 is the probability 
that the true value of 9 lies in the volume V. Given a likelihood function, the posterior 
pdf can be constructed by invoking Bayes' theorem, 

V(9\X)d9cx C(X\9) ■ n{9)d9 , (5.1) 

but the construction requires us to specify n(9), the probability density on the parameter 
space 9 prior to observing the data X. Since this prior pdf is independent of the data, 
it needs to be chosen according to one's theoretical prejudice, and is thus inherently 
subjective. 

In the often encountered situation in which no unique theoretically motivated prior 
pdf can be derived, one may wish to use one which does not favour any parameter region 
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in particular. A common choice in this case is the top-hat, or uniform, prior 
(a\AO J if #min < < # max , 

«Ue)M - | Q) otherw . se; (5.2) 

if the general order of magnitude of the parameter is known. Here, the limits 9 m m and 
# max should be chosen such that they are well beyond the parameter region of interest. 
If even the order of magnitude is unknown, one may want to choose a uniform prior in 
log 9 space instead, 

vr log (log^)dlog^=| Q; otherwige) (5.3) 

which is equivalent to a d9/9 prior in 9 space. Note that because the volume element 
d9 is in general not invariant under a parameter transformation / : 9 — >■ 9', a uniform 
prior pdf on 9 does not yield the same probabilities as a uniform prior pdf on 9' unless 
the mapping / is linear. The same is also true for the posterior probabilities, i.e., 
V(9\X)d9 ^ V(9'\X)d9' in general. 

While the posterior pdf technically contains all the necessary information for the 
interpretation of the data, the fact that it is a function in the iV-dimensional space 
of parameters makes it difficult to visualise if N > 2. Fortunately, by virtue of 
being a probability density, its dimensionality can be easily reduced by integrating out 
less interesting (nuisance) parameter directions ipi, yielding an n-dimensional marginal 
posterior pdf, 

oc Jdip 1 ...dil) m V(e 1 ,...,e n ,il) 1 ...,il) m \X), (5.4) 

which is more amenable to visual presentation if n — 1, 2, and can be used to construct 
constraints on the remaining parameters. 

A complementary approach to the marginalisation is to project the likelihood 
function C(X\0) onto the n-dimensional subspace by maximising along the nuisance 
directions, i.e, 

£prof(X|0i, -;9 n ) oc max C{X\9 U 9 n , V>i..., ip m ) . (5.5) 

W\—Wm 

Maximisation is not a Bayesian procedure, and the resulting profile likelihood cannot be 
interpreted as a probability density function. However, because £ pro f is by construction 
insensitive to our choice of priors and associated volume effects, it can be a useful 
means to assess if the inference has been significantly affected by our choice of nuisance 
parameterisation. || 

5.1. Priors 

The main parameters of interest in this work are mDM and a^ 1 . These are accompanied 
by a set of astrophysical and experiment-specific systematic nuisance parameters, as 
discussed in section 4. 

| We always normalise the (marginal) posterior pdf and profile likelihood so that max(7 : ' mar ) = 
max(£ pi . of ) = 1. 
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When it comes to specifying prior pdfs for 777dm and cr„ , we have very little guidance 
from theory without resorting to specific dark matter models. As long as the dark matter 
is cold, massive and weakly enough interacting, pretty much all combinations of values 
are a priori allowed. It thus appears reasonable to impose uniform priors on both 
log mDM and loga^ 1 . With the assumption that the dark matter particle is a WIMP, we 
can at least roughly confine our prior region. For definiteness, we take log(mDM/GeV) 
to lie in the range — > 3 and allow log(a;f /cm 2 ) to vary between —46 — > —38, as 
reported in table 1. 

Interestingly, the choice of prior boundaries on 777dm and cr^ 1 also translates directly 
to how likely we deem the direct detection experiments to actually make a positive 
detection. Consider for instance the loss of detection sensitivity for large DM masses 
(due to the large mass splitting between the DM particle and the nucleus), or for 
very light WIMPs (because of the energy threshold): the larger the prior-space in the 
{ttidm, ^j-plane, the smaller the relative fraction that the experiments will be able to 
constrain, and the smaller the subjective prior probability for them to see something. 

Our priors for the astrophysical parameters M vir and c v i r are listed in table 2. The 
ranges for the concentration parameters are inferred from simulations [74] for the NFW 
and Einasto profiles, and from fits the rotation curves of galaxies for Cored isothermal 
and Burkert profiles [75]. Note that since M vir is well-constrained by measurements (see 
table 3), the likelihood at the prior boundaries is negligible, and the inferred posterior 
pdf will be independent of our exact choice of prior boundaries for this parameter. 

5.2. Numerical implementation and construction of parameter constraints 

We employ a modified version of the public MCMC code CosmoMC [86,87], which 
uses the Metropolis-Hastings algorithm [88, 89] to sample the posterior over the full 
parameter space. The resulting chains are analysed with an adapted version of the 
accompanying package GetDist, supplemented with matlab scripts from the package 
SuperBayeS [29,90]. One- or two-dimensional marginal posterior pdfs are obtained 
from the chains by dividing the relevant parameter subspace into bins and counting the 
number of samples per bin. An x% credible interval or region containing x% of the 
total volume of P m ar is then constructed by demanding that "P mar at any point inside 
the region be larger than at any point outside. In the one-dimensional case, a credible 
interval thus constructed corresponds to the Minimal Credible Interval of [91]. Our 
profile likelihoods are also computed using CosmoMC, but with a 100-fold increase in the 
number of likelihood evaluations, so as to ensure that the tails of the distributions are 
well sampled and the true global maximum located. 

Provided the data are sufficiently constraining — that is, if the prior pdf is nearly 
constant and, under a parameter transformation f : 9 —> 9', the mapping / is almost 
linear over the parameter region where the likelihood is large — the marginal posterior 
typically exhibits very little dependence on the choice of prior. For data that can only 
provide an upper or a lower bound on a parameter (or no bound at all) however, the 
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properties of the inferred posterior and the boundaries of credible regions can vary 
significantly with the choice of prior as well as its limits m i n and 9 m3X , making an 
objective interpretation of the results rather difficult. As we shall see in the next section, 
this is in fact the case for the inference of credible regions in the {moM, cr„}-plane from 
XenonlOO, CDMSSi and CDMSGe. 

In these cases, in addition to computing credible intervals from the fractional 
volume of the marginal posterior in the {ttidm, cr^ 1 }-subspace P m ar(^DM, &^\X), we 
also construct intervals based on the volume of the marginal posterior in S'-space 
'Pm&r(S\X), where S is the expected WIMP signal, using a uniform prior on S with 
a lower boundary at zero [92]. An x% upper bound thus constructed has a well- 
defined Bayesian interpretation that the probability of S < S x is x%. The limit 
S x is then mapped onto the {mcM, a^j-plane by identifying those combinations of 
7Bdm and a^ 1 with P mar (mDM, ^I^O = 7- > ma,r(S x \X). An x% contour computed in 
this manner has the property of being independent of our choice of prior boundaries 
for mj)M and cr^ 1 . Its drawback, however, is that it has no well-defined probabilistic 
interpretation in {m DM , a^j-space.^f To distinguish these S-based credible intervals 
from the conventional ones based on the volume of P ma r(^DM, ^[AT), we label them 
with a subscript "S", e.g., 90 5 %. 

6. Results 

We present our inference results in three parts. In section 6.1 we discuss the preferred 
parameter regions in m DM and cr^ 1 for each experiment assuming the SMH (i.e., fixed 
astrophysics), after marginalising over the nuisance parameters of the experiments. In 
section 6.2 we vary in addition the WIMP velocity distribution in accordance with the 
DM density profile defined in section 3, and consider the effect of uncertainties in the 
astrophysics parameters on the inferred WIMP parameter values. Finally in section 6.3 
we entertain the possibility of a combined analysis of the DAMA and the CoGeNT data. 

6.1. Standard model halo 

DAMA Figure 1 shows our inference for the DAMA 36-bin data. The top panel 
shows the 2D marginal posterior pdf and the profile likelihood in the {ibdm, ^f}- 
subspace, where the two quenching factors gNa and q\ have been integrated and 
profiled out respectively. Both approaches single out two preferred islands of parameter 
space in {wdm,^ 1 }. Moreover, the colour coding indicates that P mar (wiDM, ) and 

% Clearly, the definition of the 5-bascd bound and its associated probabilistic interpretation are 
contingent to our choice of a uniform prior on S; Had we chosen a different prior a different set of 
limits would have resulted. Our motivation for using a uniform prior stems from the observation that, 
for Poisson statistics, a Bayesian limit on S constructed in the manner described turns out to have a 
well-defined interpretation in classical statistics, albeit a coincidental one [93]. This means the S-based 
bounds in this work can also be viewed as examples of the hybrid Bayesian/classical approach discussed 
in [94]. 
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Figure 1. Inference for DAM A assuming the SMH. Top left: 2D marginal posterior 
pdf in the {jtidm, cr^j-plane. The black solid lines enclose the 90% and the 99% 
credible regions. Top right: Profile likelihood in the {todm, c^j-plane. The black 
solid contours correspond to Axe ff = 4.6, 9.2. Bottom left: 3D marginal posterior pdf 
for {wdm,^ 1 , 9Na}, where the gNa direction is represented by the colour code. Bottom 
right: Same as bottom left, but for {tudm, j Qi}- 



£ pro f ( m DMi o"^ 1 ) coincide to an excellent degree, suggesting that the nuisance directions 
contribute no strong volume effects. For the profile likelihood, we also plot two ^-Xcs 
contours, defined via 

A Xeff( m DM,0 = -21n£ prof (m DM ,a^) , (6.1) 

where the choice of Axis = 4-6, 9.2 coincides with the classical 90% and 99% confidence 
intervals for two degrees of freedom (assuming Wilks' theorem holds). Again, we find 
remarkable agreement between these contours and the 90% and 99% credible regions 
inferred from the volume of the 2D marginal posterior. This agreement indicates that 
when the data are sufficiently informative so that the likelihood function overcomes the 
dependence on the priors, Bayesian and classical statistical methods yield very similar 
inference results. 
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Figure 2. Nuisance parameters for DAMA assuming the SMH. Left: ID marginal 
posterior pdf (black solid line) and profile likelihood (blue dashed line) for the 
quenching factor q^ a . Right: Same as left panel, but for q\. 

The bottom panel of figure 1 illustrates the correlation between {m DM , o"^ 1 } and 
the quenching factors gNa and qj. As expected, the high mass (mDM ~ 0(100) GeV) 
island is insensitive to q^, as indicated by the equal representation of gN a values in the 
island. Conversely, the low mass (mDM ~ 0(10) GeV) island shows a strong correlation 
between gN a and wdm, with higher values of gNa favouring the lower masses. The 
quenching factor for iodine shows the opposite trend: the low mass island is insensitive 
to qj, while the high mass island finds combinations of low mDM and cr^ 1 values favoured 
by large values of q\. Ultimately, however, the DAMA 36-bin data do not constrain 
either gN a or gi, as is evidenced by the fact that all values of gNa and gi allowed by their 
respective priors are represented in figure 1. The same conclusions can be drawn also 
from figure 2, which shows an essentially flat ID marginal posterior pdf (black solid line) 
and profile likelihood (dashed blue line) for gN a , while for gi one might claim a small 
preference for q\ = 0.07 — > 0.08 although it is statistically insignificant. 

CoGeNT Figure 3 shows the preferred {wdm, fff}, both in terms of the 2D marginal 
posterior pdf and the profile likelihood. As in the case of DAMA, the nuisance directions 
do not contribute strong volume effects, so that both the 90% and 99% credible 
regions inferred from the marginal posterior coincide well with the A%g ff = 4.6, 9.2 
contours on the profile likelihood surface, and single out a peak at mDM ~ 8 GeV and 
cr^ 1 ~ 10~ 40 cm 2 as the favoured region. The preferred values for the nuisance parameters 
are reported in table 4. Our analysis is compatible with all previous analyses of the 
CoGeNT data, and also with the newest data release [95], which claims detection of an 
annual modulation and where the total rate excess leads to a slightly smaller region in 
the {m DM , a^j-plane. 
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Figure 3. Inference for CoGeNT assuming the SMH. Left: 2D marginal posterior pdf 
in the {todm, cr^j-plane. The black solid lines enclose the 90% and the 99% credible 
regions. Right: Profile likelihood in the {raDM, er„}-plane. The black solid contours 
correspond to Ax^ ff = 4.6, 9.2. 



Table 4. ID marginal posterior pdf modes and 90% credible intervals for the CoGeNT 
nuisance parameters. 



Parameter 


Preferred value 


So 

C 
G n 


e.ii^ 1 kev 

4.0±t;§ cpd/kg/keV 
2.1 ± 0.5 cpd/kg/keV 



XenonlOO Our inference results for XenonlOO are shown in figure 4. Firstly, we note 
that both the 2D marginal posterior pdf and the profile likelihood form a plateau 
as m DM and cr^ 1 approach their respective lower boundaries." 1 " In this case, credible 
regions constructed from the volume of the marginal posterior in {mdm, cr^j-space 
can be strongly dependent on our choice of the and a^ 1 prior boundaries. This 

is illustrated in the left panel of figure 4 and in figure 5. In both figures the 90% 
credible region is demarcated by the black solid line, except that in figure 5 we have 
chosen a set of prior boundaries for wdm and u^ 1 (0.5 < log(mDM/GeV) < 2 and 
—45 < log(cx;f /cm 2 ) < —39) differing from the default choices of table 1. The 
discrepancy between the encompassed parameter space is clear. As an example, while 
the point {log(m D M/GeV) = 0.8, log(cr;f /cm 2 ) = —40} sits outside the 90% credible 
region in figure 4, it sits comfortably within in figure 5. 

+ Strictly speaking, the profile likelihood shown in figure 4 for XenonlOO is a quasi-profile likelihood, 
computed after the full likelihood function (4.21) has been analytically marginalised over the 
background uncertainties. 
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Figure 4. Inference for XenonlOO assuming the SMH. Left: 2D marginal posterior 
pdf in the { 772dm, ofjj-plane. The black solid line indicates the 90% bound inferred 
from the volume of the marginal posterior, while the black dashed line denotes the 
invariant 90s% contour. Right: Profile likelihood in the {ttidm, 0"„}-plane. The black 
dashed line corresponds to Ax^ ff = 2.7. 




lo 9,0( m DMK QeV ) 



Figure 5. 2D marginal posterior pdf in the {todMj <r^}-plane for XenonlOO assuming 
the SMH and an alternative set of prior boundaries for uium and cr^ 1 . The black solid 
line corresponds to the 90% bound inferred from the volume of the marginal posterior, 
while the black dashed line represents the invariant 90s% contour. 



On the other hand, the 90s% bound (black dashed line in left panel of figure 4 and in 
figure 5) is clearly independent of the boundary conditions as discussed in section 5.2, 
and the parameter region enclosed compares well with the A%g ff < 2.7 (or S < 5.2) 
region in the profile likelihood (right panel of figure 4). We will therefore use the 905% 
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Figure 6. Posterior pdf for CDMSGe assuming the SMH in the {to D m, em- 
plane. The black solid line indicates the 90% bound inferred from the volume of the 
posterior, while the black dashed lines denote the invariant 90s% and 99g% contours 
(corresponding to Ax^ ff = 3, 7.4). 



bound in the following discussion. 

Our exclusion limit on a^ 1 at high WIMP masses (wdm > 30 GeV) agrees very well 
with that provided by the XenonlOO collaboration [19]. However, at low WIMP masses, 
our bound on m DM is much less constraining compared with all previous analyses [7, 8, 
19, 25, 96]. This is clearly a consequence of the uncertainties in the scintillation efficiency 
L c ff(E) in the low recoil energy (1 < E'/keVnr < 3) region, which we have accounted for 
in this work using the nuisance parameter m. The preferred value for this parameter 
is m = 0.07 ± 0.04 (90% C.I.), which corresponds to a marginal preference for a gentler 
slope for Leff(E) at 1 < i?/keVnr < 3 with respect to the XenonlOO collaboration's 
best-fit. 

CDMSGe The posterior pdf as a function of itl-qm and a^ 1 is shown in figure 6. Since 
there are no nuisance parameters — besides the background uncertainty which we have 
already marginalised analytically in order to obtain the effective likelihood (4.7), the 
posterior pdf in figure is the full posterior pdf of the problem. It also coincides with the 
effective likelihood (4.7) because of our choice of uniform priors. A peak can be seen at 
a DM mass of 23 GeV and a cross-section of 9 x 10~ 44 cm 2 . While this is a tantalising 
hint, a detection cannot be called because the probability density is still significant at 

* We note that the exclusion limits reported by the Xenon collaboration in [19] and [66] are in fact 
ID limits on er^J for fixed values of 771dm- These limits are naturally different from our 2D limits for 
{™dm,^'}i which come from considering the joint probability distribution of jtidm and ct^ 1 . 
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much of the prior boundaries (V ~ 0.1). 

It then remains for us to set an exclusion limit in the {todm, cr„}-plane. The 90% 
contour inferred from the volume of posterior (black solid line) forms a semi-closed region 
subject strongly to our choice of prior boundaries. The invariant 90^% or Axl s < 3.0 
region (black dashed line), however, is a closed island in the {ttidm, 0"if}-plane, while 
the 99s% (Ax^ ff = 7.4) contour indicates an exclusion limit. 

Compared with the analysis in figure 3 of [24], our posterior pdf/likelihood appears 
to be more strongly peaked relative to the plateau, leading to a closed 90g% region 
while [24] finds an open one. At the same time, our peak region appears to be much 
broader than that of [24], so that their 90% contour runs right into our peak region 
where the posterior pdf/likelihood is still high (> 0.5). Fixing the number of background 
events to the mean value, i.e., setting erg = in the effective likelihood (4.7), does not 
ameliorate the discrepancy. Since reference [24] does not specify the likelihood function 
used in their analysis, we have no more handle to trace the origin of the disagreement. 

CDMSSi The analysis of the CDMSSi data is summarised in figure 7. Similar to 
CDMSGe, the posterior pdf presented in the figure is the full posterior pdf of the problem 
(barring analytic marginalisation over the background) and coincides with the effective 
likelihood (4.4). As in the case of XenonlOO, the CDMSSi data are not sufficiently 
constraining to isolate a preferred region, so that the 90% credible region inferred from 
the posterior volume (black solid line) depends on our choice of prior boundaries. On 
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lo 9lo( m DM)( GeV ) 

Figure 8. 2D credible regions for the individual experimental bounds and regions 
assuming the SMH, combined in a single plot. For DAMA (shaded) and CoGeNT 
(cyan) we show the 90% and 99% contours. The black solid line represents the 90s% 
bound for CDMSSi, and the pink dot-dash curve for XenonlOO. For CDMSGe we show 
both the 90s% and 99s% contours in blue dashed lines, while the red dotted line is 
the 90% contour for CDMSGe(LE) corresponding to Axl s = 4.6. 

the other hand, the 905% region (black dashed line, corresponding to Ax^s < 4.2 or 
S < 3.3) is independent of the m DM and cr^ 1 prior boundaries, and agrees well with the 
exclusion limit constructed by the CDMS collaboration [17]. 

SMH state of the art We summarise our results for fixed astrophysics in figure 8, in 
which we show all experimental constraints in one plot. For DAMA and CoGeNT we 
indicate the 90% and 99% credible regions, while for the exclusion limits of the other 
three experiments we show the invariant 90$% contours (also 995% f° r CDMSGe). 

We find that the parameter region favoured by DAMA is incompatible with the 
905% credible regions of XenonlOO and CDMSSi, and partially allowed by the 995% 
region of CDMSGe. In contrast, the CoGeNT preferred region is only marginally 
incompatible with these exclusion limits. Of particular interest is the compatibility 
between CoGeNT and XenonlOO. While the XenonlOO collaboration claims that their 
exclusion limit has ruled out the CoGeNT preferred region [19], we find that when 
uncertainties in the scintillation efficiency L e g(E) at low recoil energies are accounted 
for, the CoGeNT and the XenonlOO data can find some common ground. 
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Between CoGeNT and DAMA we find that their 99% credible regions do not 
overlap, despite marginalisation over the quenching factors q^a and qi for DAMA. This 
is a consequence of our choice of prior boundaries for q^a. (0.2 — > 0.4), especially in 
view of [6, 7], where it has been suggested that in order to make DAMA and CoGeNT 
compatible large quenching factors for sodium (e.g., gNa — 0.6) and for germanium 
should be considered. Allowing up to 10% of channelling for DAMA could also improve 
the agreement between the two experiments, by shifting the DAMA low mass region 
downwards in the {wdm, a^j-plane [69]. 

Lastly, we also show the exclusion bound derived from CDMSGe(LE) (red dotted 
line in figure 8). Since the likelihood function (4.11) for this case is a multivariate 
gaussian, we can infer an invariant 90% exclusion bound similar to the 90s% bound 
by demanding that Ax^ ff < 4.6. At low masses this 90% bound turns out to be very 
close to the CDMSSi exclusion limit, so that the DAMA preferred region falls outside 
the credible region, while the CoGeNT region falls mostly within. The main difference 
between this low energy analysis and the standard CDMSGe is that the former does not 
find any closed region at small WIMP masses. Compared with other experiments, for 
masses larger than 10 GeV the XenonlOO bound is more constraining. We therefore do 
not consider the CDMSGe(LE) exclusion limit any further. 

6.2. Variable WIMP velocity distribution and astrophysics 

Figure 9 shows the effects of astrophysical uncertainties on the inferred {m D M, o"^ 1 } 
parameter space, for each of the four DM density profiles considered (see section 3.1). 
The corresponding preferred values for the local dark matter density, the circular and 
the escape velocities are reported in table 5. 

Firstly, we note that all four DM profiles give very similar inference results on the 
{mDM, cr„}-plane. This means that the exact shape of the DM halo density profile — at 
least within the class of spherically symmetric, smooth profiles — does not yet play a 
role in direct DM searches. This conclusion is further supported by the inferred local 
DM density, circular and escape velocities presented in table 5. The preferred values for 
these quantities differ from profile to profile, with the Cored isothermal halo in particular 
favouring the very high end of the observationally allowed escape velocities (see table 2). 
However, once the DM halo profile has been fixed, we see that the preferred values for t> , 
v esc and p Q and their associated uncertainties are virtually independent of the additional 
constraints from the DM experiments. In other words, direct DM searches are not at 
the moment contributing towards constraining the astrophysics of the problem. 

Secondly, we note that allowing for uncertainties in the astrophysics significantly 
expands the closed regions of DAMA, CDMSGe and CoGeNT, while the exclusion limits 
tend to shift a little to the right. For all four profiles, the preferred regions of DAMA 
and CoGeNT now appear to marginally overlap: for the NFW, Einasto and Burkert 
profiles we see an overlap between the 90% credible region of DAMA with the 99% 
region of CoGeNT and vice versa, while for the Cored isothermal the agreement is a 
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Figure 9. Same as figure 8, but with variable astrophysics, assuming the Cored 
isothermal (top left), NFW (top right), Einasto (bottom left), and Burkert (bottom 
right) profiles. 
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Figure 10. 3D marginal posterior pdf for DAMA and CoGeNT for {todm,^ 1 } and 
the circular velocity vo (left), the escape velocity v csc (centre), and the local DM density 
Pq (right), assuming the NFW profile. The third parameter direction is represented 
by the colour code. 
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Table 5. ID posterior pdf modes and 90% credible intervals for the circular velocity 
vo, escape velocity z; eS c, and the local DM density p for DM density profiles considered 
in this work. 
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little worse. One may be tempted to claim some degree of agreement between DAMA 
and CoGeNT based on this partial overlap. However, before we do so, it is important 
that we also examine the degree of overlap between the preferred regions in the other 
parameter directions. 

Figure 10 shows the 3D marginal posterior pdf for {mDM, c^ 1 } and a third parameter 
direction t>o, f eS c and p . Here, we see that while it is not impossible to find a value of 
f csc that satisfies both DAMA and CoGeNT simultaneously, there is a clear trend that 
combinations of larger {m^M, } values tend to prefer higher values of Vq (and similarly 
for p ). This indicates that although DAMA and CoGeNT appear to overlap in the 
{rooM, cr^ I }-plane, there is in fact very little overlap between them in the v direction 
(and naturally also in the p direction which enters into the differential recoil rates as 
a common normalisation factor). 
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Figure 11. Inference from the combined DAM A and CoGeNT fit assuming the 
Burkert halo profile. Left 2D marginal posterior pdf in the {ttidm, <7„}-plane. The 
black solid lines correspond to the 90% and 99% bound inferred from the volume of 
the marginal posterior. Right: ID marginal posterior pdf (black solid line) and profile 
likelihood (blue dashed line) for the quenching factor q^ a . 

6.3. Combined fit? 

Despite apparent difficulties to reconcile the DAMA and the CoGeNT preferred 
regions within the boundaries of our nuisance and astrophysics models, let us for a 
moment entertain the possibility of a combined fit. Figure 11 shows the preferred region 
in the {ttidm, cr„ }-plane from a combined fit of DAMA and CoGeNT assuming the 
Burkert profile (marginalised over all nuisance and astrophysics parameters as usual). 
The corresponding ID marginal credible intervals for the nuisance and astrophysics 
parameter are displayed in table 6. 

The best-fit point of the combined fit corresponds to a mass of 9.2 GeV and a 
cross-section of 1.26 x 10 -40 cm 2 . However, this fit comes at the expense of a significant 
shift in the circular velocity: Vq = 176^^ 3 km s _1 (90% C.I.) from the combined fit, 
versus vq = 214^2? km s -1 from fitting either DAMA or CoGeNT alone. The preferred 
local DM density p Q and escape velocity v esc also suffer a downward shift respectively, 
although not a significant one in either case. For the nuisance parameters, we find that 
for CoGeNT, the normalisation for the exponentially decaying background C has come 
down a little, and the decay rate Sq is significantly more constrained. The radiative 
peaks, on the other hand, are a well-defined feature, and consequently the preferred 
value for their height G n has not been affected by the combined fit. 

Most interestingly, we find that the DAMA sodium quenching factor gNa, which 
was previously an unconstrained quantity (see figure 2), now shows a preference for high 
values (right panel of figure 11). In particular, the ID profile likelihood (blue dashed 
line) has no local maximum, and hits its highest point right at the prior boundary 
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Table 6. ID marginal posterior pdf modes and 90% credible intervals for the 
astrophysical and the nuisance parameters from the combined DAMA and CoGeNT 
fit assuming the Burkert profile. 



Parameter Preferred value 



^0 


176tf km s" 1 


^esc 


533tf km s" 1 


P© 


0.3 + ^ 9 GeV cm" 3 


<?Na 


o.3818:8§ 


So 


5 ± 1.2 keV 


C 


2.81?;? cpd/kg/keV 


G n 


2.2 ± 0.4 cpd/kg/keV 



<?Na = 0.4. This suggests that if we had allowed for a wider prior range for g Na , an 
even higher value might have been preferred. This result is consistent with previous 
suggestions that a higher value for g Na could improve the compatibility of DAMA and 
CoGeNT (see section 6.4). 

To assess the quality of the fit for the best-fit point singled out by the combined 
run, we look at the spectral shape of the expected signal in both detectors. In figure 12 
we show on the left the averaged modulated amplitude of DAMA, and on the right 
the number of counts per bin versus the recoil energy for CoGeNT. Superimposed here 
are the predictions for the best-fit point (dashed lines). Clearly, the "best-fit" point 
is actually a bad fit for both experiments. For the best-fit DM mass, cross-section 
and nuisance parameters, a better fit in DAMA would be obtained by increasing both 
t>o, which would shift the spectral curve to the left, and p , which would result in a 
global enhancement of the signal. For CoGeNT the trend is the opposite: a better fit 
is obtained by decreasing p and increasing t>o, as demonstrated in figure 10. For both 
detectors the signal is rather insensitive to the value of v esc . 

6.4. A larger g Na 

As suggested in the previous section, allowing for a larger sodium quenching factor g^a 
for DAMA may improve the combined DAMA/CoGeNT fit. We explore this possibility 
here by raising the upper limit of our prior range on q^^ from 0.4 to 0.6, and recomputing 
the preferred regions for combined DAMA/CoGeNT assuming a Burkert profile. 

The results are shown in figure 13. On the left panel, the preferred region in 
the {wdm, <7^ 1 }-plane is similar to that inferred using our standard prior on g Na (see 
figure 11), with the best-fit point now corresponding to a mass of 7.38 GeV and a cross- 
section of 9.64 x 10~ 41 cm 2 . On the right panel, we see that both the ID marginal 
posterior pdf (black solid) and profile likelihood (dashed blue) for g Na rise sharply 




Figure 12. Left: The expected signal for DAMA using the best-fit point of our 
combined DAMA/CoGeNT fit. Right: The expected signal for CoGeNT using the 
best-fit point of our combined fit. 



between = 0.5 and 0.6, hitting their highest points at the edge of the prior boundary. 
This confirms the trend that the combined DAMA and CoGeNT data prefer a large q^. 

Interestingly, the preferred values for the astrophysical parameters from the 
extended combined fit are now more closely in line with those from fits to the individual 
experiments alone (see table 7). The most significant change can be seen for the 
circular velocity, which now has the preferred value vq = 201^7 km s _1 (90% C.I.), 
in contrast to (a) v = 176^f 3 km s -1 from the standard combined fit in section 6.3, 
and (b) v = 2U±H km s" 1 and v = 216±| km s" 1 from DAMA and CoGeNT alone 
respectively. 

However, despite this shift in the astrophysical parameters, the extended combined 
fit offers only a marginal improvement over the standard combined fit. This can be 
seen in figure 14, where we show the spectral shapes of the expected signals for the 
individual experiments corresponding to the best-fit point of the extended combined fit. 
Comparing with figure 12, we see that the fit to CoGeNT now shows better agreement 
to the data at low energies, while for DAMA the higher value for gNa now leads to a 
lower peak in the spectrum, which is further suppressed by the smaller value of cr^ 1 and 
the lighter dark matter mass. 

7. Conclusions 

The present status of the direct detection of dark matter is somewhat ambiguous. On the 
one hand, there have been claims of detection of a low-mass WIMP signal from DAMA 
and CoGeNT. On the other hand, the XenonlOO, CDMS, and CDMS-II experiments 
have only been able to provide exclusion limits. The interpretation of these experiments 
in terms of a dark matter signal is complicated by the presence of backgrounds, and the 
need to model experiment-specific systematic effects, such as the quenching factors for 
DAMA or the scintillation efficiency for the Xenon detector. In addition, one requires 
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Figure 13. Same as figure 11, but for an extended prior range for the DAMA sodium 
quenching factor QNa (up to gNa = 0.6). 



Table 7. Same as table 6, but for an extended prior range for <?Na (up to (?Na = 0.6). 



Parameter Preferred value 

v 20lip km s" 1 

v esc 54]+^ km s" 1 

p Q 0.36i{j:og GeV cm" 3 

Qnsl 0.59lg;o4 

S 9.4 ± 1.8 keV 

C 3.1±?;1 cpd/kg/keV 

G n 2.2 ± 0.4 cpd/kg/keV 



the input of astrophysical quantities, which enter into the theoretical expressions for 
the total or modulated DM rates. All these effects are to some extent subject to 
uncertainties, which need to be propagated to the inferred dark matter parameters. 

This multi-parameter inference problem can be addressed in a simple and consistent 
way using Bayesian statistical methods. In the present work, we apply these methods 
to a selection of current direct dark matter searches to infer the mass and cross- 
section of WIMP dark matter in the simplest scenario of spin-independent elastic WIMP 
scattering. 

We initially ignore the astrophysical uncertainties and focus on the effects of 
experimental nuisance parameters and background uncertainties. Our main result is 
that the XenonlOO exclusion bound is significantly weakened once the uncertainty on 
the scintillation efficiency is taken into account. As a consequence, we find that the 
CoGeNT preferred region in the {m DM , a^j-plane is quite compatible with XenonlOO, 
and there is even a marginal consistency (at 90s% credibility) with the DAMA preferred 
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Figure 14. Same as figure 12, but for the best-fit point of the extended combined fit, 
with mrjM = 7.38 GeV, cr„ = 9.64 x 10~ 41 cm 2 , and other parameter values presented 
in table 7. 



region. We expect that this conclusion also holds for XenonlO and the latest CoGeNT 
data which were obtained after two years of data taking [97]. We also remark that 
after marginalising over the background, the standard CDMSGe analysis yields a closed 
90s%-credible contour, although no closed regions remain in the low energy re-analysis 
CDMSGe(LE). 

We then repeat the analysis procedure including astrophysical uncertainties, 
considering besides the standard model halo three other spherically symmetric dark 
matter halo models with isotropic velocity distributions, whose density profiles are 
motivated by iV-body simulations. We find that the inferred values of the astrophysical 
parameters are independent of the direct detection experiment data, indicating that 
their values depend only on the chosen DM density profile. With the exception of 
the isothermal halo, which prefers significantly higher escape velocities, the different 
halo parameterisations lead to similar values of the astrophysical parameters. Not 
unexpectedly, including the astrophysical uncertainties further reduces the constraining 
power of the data in the {itidm, cr^j-plane. At first glance, this seems to improve 
the compatibility between DAMA and CoGeNT. However, this impression is somewhat 
misleading, since in some of the marginalised directions (vq and p & ) the disagreement 
remains — tweaking astrophysics parameters alone cannot reconcile the two results. If we 
demand compatibility between these experiments, then the inference process naturally 
concludes that a high value for the sodium quenching factor for DAMA is preferred. 

It will be interesting to apply the analysis framework presented in this paper to 
more complex models of the dark matter halo, like asymmetric velocity distributions or 
the presence of streams in the Galactic halo. Additionally, an application to alternative 
scenarios for the particle physics interactions can be envisaged, such as inelastic DM [4, 
95] or more exotic scenarios, as for instance discussed in [4,24,98]. 
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Appendix A. Dark matter density profile in terms of M vir and c cir 

The two-parameter DM density profiles defined in section 3.1 in terms of p s and r s can 
be expressed in terms of the halo's virial mass M vir and concentration parameter c v ; r . 
Firstly, the parameter r s can be parameterised as 

r s (M vir ,c vir ) = rvir(Mvir) , (A.l) 

Cvir 

where the virial radius r vir defines a spherical region in which the average DM density 
is S c = 200 times the critical density p crit . The mass enclosed in this region is called the 
virial mass, 

M vir = An f ™ dr r 2 p DM (r) = ^7rr vir 5 c p crit . (A.2) 
Jo 3 

Using this relation, we can solve for p s once a profile has been specified. We give the 
solutions for the four halo profiles considered in this work: 

(i) Cored isothermal : 



(ii) NFW: 



(iii) Einasto: 



^-^^ (A.3) 



r 3 

3 ln(l + c vir ) - c vir /(l + c vir ) 



( v _ a c Pcrit<4[2 a exp(|)ft^ j 1 

where T(a) and F(a, b) are the gamma and the incomplete gamma functions, 
respectively. 

(iv) Burkert: 

/ \ 4<5 c p cr it c vir . . 

Psl ° virJ= 3 2 ln(l + c vir ) + ln(l + 4 r ) - 2 tan^c^) ' 1 ' 
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